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Abstract 

A nonlinear nonlocal cochlear model of the transmission line type is studied in order to 
capture the multitone interactions and resulting tonal suppression effects. The model can 
serve as a module for voice signal processing, it is a one dimensional (in space) damped 
dispersive nonlinear PDE based on mechanics and phenomenology of hearing. It describes 
the motion of basilar membrane (BM) in the cochlea driven by input pressure waves. 
Both elastic damping and selective longitudinal fluid damping are present. The former 
is nonlinear and nonlocal in BM displacement, and plays a key role in capturing tonal 
interactions. The latter is active only near the exit boundary (helicotrema) , and is built in 
to damp out the remaining long waves. The initial boundary value problem is numerically 
solved with a semi-implicit second order finite difference method. Solutions reach a multi- 
frequency quasi-steady state. Numerical results are shown on two tone suppression from 
both high-frequency and low-frequency sides, consistent with known behavior of two tone 
suppression. Suppression effects among three tones are demonstrated by showing how 
the response magnitudes of the fixed two tones are reduced as we vary the third tone in 
frequency and amplitude. We observe qualitative agreement of our model solutions with 
existing cat auditory neural data. The model is thus simple and efficient as a processing 
tool for voice signals. 
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1 Introduction 



Voice signal processing algorithms have received increasing attention in recent years for im- 
proving the design of hearing devices and for evaluating acoustic theories of peripheral auditory 
systems, see Meddis et al |23 among others. A fundamental issue a computational method 
needs to resolve is the nonlinear interaction of acoustic waves of different frequencies and the 
resulting tonal suppression effects. A nonlinear filter bank approach is developed in based 
on knowledge of auditory responses. Nonlinearities are known to originate in the cochlea and are 
further modified in higher level auditory pathways. The cochlear mechanics has first principle 
descriptions, and so partial differential equations (PDEs) form a natural mathematical frame- 
work to initiate computation. However, in vivo cochlear dynamics is not a pure mechanical 
problem, and neural couplings are present to modify responses. To incorporate both of these 
aspects, we aim to develop a first-principle based PDE approach to voice signal processing, 
where the neural aspect is introduced in the model phenomenologically. 

Cochlear modeling has had a long history, and various linear models have been studied at 



length by analytical and numerical methods, see Keller and Neu Wm, Leveque, Peskin and Lax 



lq] , Lighthill [fLql , and references therein. It has been realized that nonlinearity is essential for 
multitone interactions, see [11, [15], 0, |] among others. The nonlinearity could be incorporated 
through micro-mechanics of cochlea, such as coupling of basilar membrane (BM) to inner hair 



cells fl3| , |T9|| . Or nonlinearity could be introduced phenomenologically based on spreading 
of electrical and neural activities between hair cells at different BM locations suggested by 
experimental data, see Jau and Geisler [1^], Deng ||. The latter treatment turned out to be 
quite efficient and will be adopted in this paper. Their model with nonlinear and nonlocal BM 
damping will be our starting point. 

Multitone interaction requires one to perform numerical simulation and analysis in the time 
domain as the resulting time dependent problem is strictly speaking irreducible to a steady state 
problem. The commonly used cochlear models, including those of |12] and ||, are intrinsically 
dispersive. In particular, long waves tend to propagate with little decay from entrance point 
(stapes) to the exit point (helicotrama) . Such an issue can be avoided in linear models because a 
reduction to the steady state (or a frequency domain calculation using time harmonic solutions) 
is available by factoring out the time dependence of solutions. The dispersive phenomenon 
prompted us to incorporate a selective fluid longitudinal damping term in the model of Jau and 
Geisler |12j], Deng ||, operating near the exit point. The role of such a term is to suck out the 
long waves accumulating near the exit. Selective positive or negative damping has been a novel 
way to filter images in PDE method of image processing, see Osher and Rudin for the first 
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work in this direction. Our treatment here is similar in spirit though on a different type of PDE. 

We then present a second order semi-implicit finite difference method of our model. The 
method relaxes the stability time step constraint of explicit methods without going into the 
complexity of fully implicit schemes. As our model is nonlocal and nonlinear, semi-implicit 
scheme is a reasonable option for computation. Numerical examples showed that the method 
is efficient, and that the selective damping term indeed removes the long waves at the exit 
boundary point for a time domain calculation of a single tone. The examples also demonstrated 
that the analysis of dispersive waves and slow decay of long waves is robust, and the phenomenon 
persists in the nonlinear nonlocal regime of the model. 

Subsequently, we show computational results on isodisplacement curves for five characteristic 
frequencies consistent with those of Neely [^] for BM displacement of 1 nm. In case of two tone 
interactions, our model gives both low and high side suppression, in qualitative agreement with 
similar studies of Geisler ||. For understanding three tone interactions, we fix the frequencies 
and amplitudes of two tones, and vary the third tone both in frequency and amplitude in 
a manner similar to auditory neural experiment of Deng and Geisler Again, our model 
showed suppressing effect on the two tones by the third tone, in qualitative agreement with 
experimental data. To our best knowledge, the current paper is the first on time domain 
simulation of interactions of three tones using first principle based PDEs. The stage is set for 
further computation, analysis, and model development on more complicated tones. 

The rest of the paper is organized as follows. In section 2, we present our model and its 
background, also analysis of dispersive properties of solutions. In section 3, we outline our 
semi-implicit discretization of the model. In section 4, we give model parameters, and show 
numerical results on dispersion of model solutions, isodisplacement curves, two and three tone 
interactions in comparison with existing data. The conclusions are in section 5. 



2 Nonlinear Nonlocal Cochlear Model 
2.1 Background and Model Setup 

The cochlear modeling has had a long history driven by advancement of cochlear experiments, 
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25, ^5], |T], [| among others. A brief derivation of the one dimensional (1-D) cochlear model 
of transmission line (long wave) type, based on fluid mechanics and elasticity, can be found 
in Sondhi ||26|| , see also references therein for earlier contributions, and de Boer || for later 



development and higher dimensional models. A general form of the 1-D model can be expressed 
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as: 

p xx - Nu a = e s u t , x e (0, L), (2.1) 
p = mu t t + r(x, u)u t + s(x)u, (2.2) 

where p is the fluid pressure difference across the basilar membrane (BM), u the BM displace- 
ment, L the longitudinal length of BM; stapes is at x = 0, and helicotrema is at x = L; N a 
constant depending on fluid density and cochlear channel size; e s depends on the damping of 
longitudinal fluid motion [20|; m, r, s are the mass, damping, and stiffness of BM per unit area, 



with m a constant, s a known function of x, and r a function(al) of x, u. 

Choosing a functional form for r is a key step in model design. Many choices in the literature 
assume that r is a local nonlinear function, [|ll|, [7|, |29|, [H| 0. The nonlinearity is necessary to 
model two tone nonlinear interaction (e.g. two tone suppression effects). However, as pointed 



out by Kim [15|, local form of r is not sufficient to account for two tone suppression when sup- 



pressors are below the excitors in frequency (low side), unless a "second filter" is incorporated, 
see e.g [11]. To achieve the second filter effect on more physical grounds, micromechanical ap- 
proaches have been developed to derive damping r from motion of outer hair cells and other 
physiologically relevant parts of cochlear, JT^, |29], ^ [Tj| among others. Though this is a rea- 



sonable approach, one faces the daunting complexity of the cochlear micromechanics, see [19] 
for recent progress along this line. The alternative approach, first proposed in this context by 
Jau and Geisler is phenomenological in nature. It hypothesized that there is longitudinal 



coupling in the cochlear partition itself, and damping r is a convolution integral of u with an ex- 
ponential weight function on BM. Similar idea appeared also in Lyon ETJ in coupled automatic 



gain control framework. Possible physiological mechanisms underlying the spatially coupled 
nonlocal damping were identified in Deng ||, and the resulting model was then used to process 
acoustic signals. 

We shall adopt the latter nonlocal damping approach for choosing the function r. The BM 
damping coefficient is specified as in Jau and Geisler fl2| , Deng 0: 



r(x,t,u) = r (x) + 7 / \u(x , t)\ exp{— \x — x l/A} dx', (2.3) 

Jo 

where tq(x) is the passive BM damping; 7 and A are positive constants. The size of r controls 
how sharp the BM traveling wave envelope would be. 

We shall take e s as a function of x so that e s = e s (x) > is a smooth function and is 
supported near the vicinity of x = L, for instance: 

es{x) = l + exp{(3(x s -x)y (2 ' 4) 
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where x s is a constant slightly smaller than L; e > 1, (3 > 1, are two constants, with /3 suitably 
large. Taking e s as a spatially dependent function can be considered as a way to selectively 
introduce damping so that possible low frequency waves near x = L are damped out and there 
is minimal wave accumulation (or reflection) close to the helicotrema (at x = L). This turns 
out to be an essential stabilizing mechanism for our time domain numerical computation. We 
will come back to this point later. 

The coefficients m, n, and s(x) are standard. The function s(x) is based on the data of 
Liberman |T7| ]: 

s(x) = 4vr 2 m (0.456 exp( 4.83 ( 1 - x/3.5 )) - 0.45 ) 2 . 
The physical boundary and initial conditions are: 

Px {0, t) = T M p T (t), p{L, t) = 0, (2.5) 
u(x, 0) = ut(x, 0) = 0, (2.6) 

where pr{t) is the input sound pressure at the eardrum; and T M is a bounded linear operator 
on the space of bounded continuous functions, with output depending on the frequency content 
of pr(t). If pt = J2j=i Aj exp{iujjt} + c.c, a multitone input, c.c denoting complex conjugate, 
Jm a positive integer, then TmPt^) = J2j=i Bj exp{iujjt} + c.c, where Bj = aM{uJj)Aj, c.c for 
complex conjugate, «m (') a scaling function built from the filtering characteristics of the middle 



ear. Established data are available in Guinan and Peake [10|. The model setup is now complete 



2.2 Model Properties and Analysis 

We briefly go over some special regimes and solutions of the model system to illustrate the basic 
mathematical properties and issues. If m — 0, r = 0, e = 0, (|2.1| )-( fE2"l) reduces to a second order 
wave equation: nu u — (s(x)u) xx = 0. When e is turned on, and j3 = 0, we have a damped second 
order wave equation driven by boundary condition at x = 0. Now suppose that all coefficients 
are constants, then (|2.1|) - (|2.2| ) considered on the entire line admits planar wave solutions of the 
form: 

u = uoexp{i(kx — cut)}, p = poexp{i(kx — out)}, 

where k is the spatial wave number, u> the temporal frequency; u and p are respectively the 
amplitude of displacement and pressure. Upon substitution, it follows that: 

Po = {—muu 2 — iruu + s)uq, 

(k 2 m + N)ou 2 + i(k 2 r + e)uu - sk 2 = 0. (2.7) 
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If the damping coefficients r = e = 0, then: 



UJ = ± (k 2 m + Ny/^ (2i 



showing that the system is dispersive, see ||30|| , i.e, to = u{k) and u"(k) ^ 0. Waves of different 
wave length (27c/k) travel at different velocities. The dispersion relation ( |2.8[ ) gives the group 
velocity: 

u'{k) = ±sl /2 N(mk 2 + N)~ 3/2 , (2.9) 

which decays to zero as k — > oo. In other words, short waves (k large) do not disperse as fast 
as long waves (k small) . 

When damping coefficients r, e > {(3 = 0), we have a damped dispersive system, Q2.8| ) is 
modified into: 



-i{k 2 r + e) ± J As k 2 (k 2 m + N) - (k 2 r + e) 2 

W = 2(Pm + iV) ' (2 ' 10) 

with Im{uj} = — (k 2 r + e s )/2(k 2 m + N) < 0, if the discriminant is nonnegative; otherwise: 



(k 2 r + e s ) ± J-Ask 2 (k 2 m + N) + (k 2 r + e s ) 2 
Im{uj} = 1 — — — — < 0. 

In either case, we see decay of planar waves. The decay rate shows that if e s = 0, the 

decay of long waves (\k\ 1) is very slow. Hence the decay of long waves relies on e s . 

We remark that the cochlear model in || || also contains a so called stiffness coupling term 



in equation Q2.2Q , so the elastic response becomes: p = mu t t + r(x } t, u)u t + s(x)u — K(x)u xx , for 
some positive function K(x), rendering the highest spatial derivative order four in the model. 
Nevertheless, it is easy to check that the additional term does not change the dispersive nature 
of model solutions, neither the long wave dispersion and slow decay properties. 

As we will see, the long wave dispersion and slow decay properties persist when coefficients 
become variables, and could cause problems for numerics. However, the generic dispersive effect 
of the time dependent solutions is absent in the special case when the temporal dependence of 
solution can be factored out. This is the case when nonlinearity is absent, i.e. r = r(x). Let the 
input Pr(t) be a single tone, or in complex variables Pr{t) = Aexp{iut} + c.c. Writing solution 
as u = U(x) exp{iut} + c.c, p = P{x) exp{iut} + c.c, we derive the following boundary value 
problem: 

P{x) = (-mu 2 + ir{x)uj + s(x))U, (2.11) 
P xx + (Nlu 2 - e(x)iuj)U = 0, (2.12) 
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subject to the boundary conditions: 



P x (0) = T M A, P(L) = 0. 



(2.13) 



This is the way frequency domain calculation is done, turning a time dependent problem into a 
steady state problem on wave amplitudes, avoiding transients and the dispersion effect. How- 
ever, when nonlinearity is present, it is in general impossible to factor out the time dependence 
exactly, especially when the input signal consists of multiple tones. To do an accurate computa- 



boundary value problem in time until transients die out and solutions approach a limiting time 
dependent state (quasiperiodic in general for multiple tones). 



numerical issues. Solutions of cochlear model with tonal inputs behave like traveling waves 
during early time, and eventually settle down to time periodic states of different periods at 
various BM locations according to the resonance information (frequency to place mapping) 
encoded in the function s(x). A time domain computation will provide the stabilized limiting 
multi-frequency oscillating states after a long enough time integration. An ideal numerical 
method would be one that can reach the limiting states with accuracy, low cost and speed. 

Let us first briefly review some existing methods. For a cochlear model with locally nonlinear 
damping function r (van der Pol type), a method of line discretization is implemented in |7|. 
Discretization in space is a Galerkin finite element method (of second order, with piecewise 
linear basis), and time stepping is explicit variable step 4th order Runge-Kutta, found to be 
superior to earlier explicit time stepping schemes such as Heun's and modified Sielecki methods. 
A fully discrete explicit finite difference discretization, first order in time and second order in 
space, is carried out in 0, 0]. The discretization can be viewed as central differencing in space 
and forward Euler in time. The explicit nature of these methods put a severe stability constraint 
on the time step, requiring long computation time to reach stable quasi-periodic states that we 
are most interested in. Also it is known that one step explicit methods such as Runge-Kutta, 
Heun and Euler, are prone to accumulation of truncation errors for approximating dynamic 
objects like limit cycles, see |J, chapter 2. 

Though implicit methods typically relax the stability constraints on time step and are better 
for computing steady states, they are also likely to be slow due to Newton iteration at each time 
step, and the non-banded complicated Jacobian matrix as a result of the nonlocal damping term. 



tion directly on system ( |2.1|) - (|2.6|) , one has to proceed in the time domain, or solving the initial 



3 Numerical Method 



In this section, we discuss our numerical method 




related 
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When r is only Lipschitz continuous in it, lack of smoothness is another concern for the condition 
of the associated Jacobian matrix. This motivates us to look into semi-implicit methods that 
allow larger time steps and better stability properties than explicit methods, while being cheaper 
and faster than fully implicit schemes. 

Let us introduce our 2nd order accurate semi-implicit discretization. We denote by itf 
(p™) the numerical approximation of solution u(jh,nk) (p(jh,nk)), where h the grid size for 
x G [0, L], k the time step. Let 6± h be the forward/backward finite differencing operator in x, 
8±k the forward/backward finite differencing operator in t. Central first differencing operator in 
space is 5 ,h = + '% ~' h , and central second differencing operator in space is 5\ = (S+^—S-.^/h; 
with similar notations <5o,fc> &\ f° r central differencing operator in time. The semi-implicit method 
for system (|2.1[ - |2.6| ) is: 

NSlu] = + 2p] + p]- 1 ) - eA*«" , (3-1) 

2< +1 - 5< + 4<~ x - u r r 2 
Vf = 3 —^ 3 — + rj(u t )] +1 + Sj u] +l , (3.2) 

where 1 < j < J, n > 2, e 3 — e s (xj), and: 



r™ = [r 0tj + -fh(-\u r l +1 \e'^- Xll/X + -\u n j +1 \e' lx '- Xjl/X + ^ |< +1 |e" |Xj " Xil/A )], (3.3) 
2 2 l=2 



where {u t )!j +l is approximated at second order by: 

(u t )^ +1 = (u t )]- 1 + (2k)(u tt )]- 1 + 0(k 2 ) 

9k 



and: 



2k 

5 Q , k u n f l + j^iSlp]- 1 - e^ k u n ~ l ) + 0(k 2 ), (3.4) 



K +1 I = K + k( Ut )? + o(k 2 )\ 

= \uf + (Au-' 1 - ur 2 - 3<)/(-2)| + 0(k 2 ) 

= |(4«r 1 -<- 2 -5<)/(-2)| + 0(A ; 2 ). (3.5) 



In ( |3.4|) , we have used the equation Q2.1|) once to lower time derivative of it by one order. In 
fl3.2|) , ( p.3|) and ( p.5|) , we have also used one sided second order differencing to approximate u a , 
u t , and trapezoidal rule to approximate the spatial integral. The variables at t = (n + l)k are 
all linear. If follows from (B.2I) that: 



un+1 = k 2 pf l - m(-5 M ? + Aur 1 ~ ur 2 ) + *x?(4tff - u^)/2 
3 2m + 3Ax? /2 + k 2 Sj ' 1 ' ' 
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Substituting ( |3.6| ) into ( |3.1| ) to eliminate u n+1 , we get the following linear system of equations 
onp n+i ( A = k / h y 



c b 1kr n \ 2 \ 2 n n+1 \ 2 r) n+l 

f )*>/(*» + ^ + *».,) - ^-bf 1 + ^ + 



ra— 1 „ n—2\ | / /i „ n „ n-\\ 

— 1 1 1 \ — ; i / / . — i— i / 

iv(-2<+ M r i ) + (iv+^ r )- 



Gj k -m(-5u] + Au]- 1 - ) + %t^ u l ~ U D 



(2m + Y + k 2 Sj ) 



--[2(p] +1 - 2p] + p]_ x ) + p]~l - 2p n f l + p]ll} - ^u]-\ (3.7) 

which is solved by inverting a diagonally dominant tridiagonal matrix. Once p n+l is computed, 
u n+1 is updated from ( |3.6|) . 

The first two time steps are initiated as follows. Initial condition, ( |2.1|) , and ( |2.2|) give: 

N 

Pxx - ~P = 0, p a = /(0), p(L, 0) = 0, 
m 

whose solution is: 

p(x, 0) = /(0) ? =^(exp{(x - 2L) v /iVy^} - expj-V^x}). (3.8) 

At £ = k, we found 

1.2 

k) = k 2 u tt /2 + 0(fc 3 ) = —p(x, 0) + 0(A; 3 ); 

and similarly: 

u t (x, k) = kp(x, 0)/m + k 2 u m {x, 0)/2 + 0(A; 3 ). 

To find p(x, k), denote q(x) = Pt(x, 0) = muut{%, 0) + r(x, 0)u tt (x, 0). It follows from ( |2.1| ) that 
q xx - Nu m (x, 0) = e(x)u tt {x, 0), implying 

N N 

Qxx q= -e(x)u tt (x, 0) r(x, 0)u tt (x, 0), 

m m 

or a uniquely solvable two point boundary value problem on q: 

qxx -—q= -(Nr(x, 0)/m 2 - e(x))p(x, 0), (3.9) 
m 

with boundary condition: q x (0) = f'(Q), q{L) = 0. 
It follows then from the definition of q that: 

u m (x,0) = q(x)/m — r(x, 0)p(x, 0)/m 2 , 

k 2 

Ut(x,k) = kp(x, 0)/m + — (q(x)/m — r(x, 0)p(x, 0)/m 2 ). (3.10) 
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Equations (|2.1|) - ([2.2|) at t — k implies the uniquely solvable two point boundary value problem 
on p(x, k): 

N 

Pxx V — ( e ( x ) ~ Nr(x,u))ut — Ns(x)u/m, (3-H) 

with boundary condition: p x (0) = f(k), p(L) = 0. Both ( ^.9[ ) and ( |3.11| ) are solved with a 
standard second order discretization. We then will have computed (u,p)(x,k). Similarly, with 
k replaced by 2k, we compute (u,p)(x,2k). 

The left boundary condition p x (x, 0) = f(t) is discretized as: p n 2 +l -Pv +1 = 2hf((n+l)k), 
with second order accuracy. The right boundary condition PjW = is exact. The semi- implicit 
method is now completely defined, and it reduces to an unconditionally stable implicit method 
for standard second order wave equation when m = 0, s is constant, and damping is absent, 
27|, chapter 8. 



4 Numerical Results and Multitone Interaction 

In this section, we present numerical results based on our model and numerical methods dis- 
cussed in the last two sections. First we list all parameters of our cochlear model and those 



of the associated middle ear model based on Guinan and Peake [IU|. We then show computed 
single tones and demonstrate the dispersion effects of long waves and the role of selective lon- 
gitudinal fluid damping in our model. As a test of model tuning property, we also show the 
computed isodisplacement curves (at characteristic frequencies 500 Hz, 2 kHz, 4 kHz, 6 kHz, 10 
kHz). Subsequently, we give numerical examples of both the high side and low side suppression 
of two tone interaction in model solutions consistent with Figures 10.2, 10.4, 10.7 in |§. The 
asymmetry of high and low side suppression is captured by the model. We then show numerical 
simulation of three tone interaction, presenting results on how the amplitudes of the two of 
the three tones change as we vary the amplitude and frequency of the third tone, qualitatively 
consistent with experimental data of Deng and Geisler 0] on responses of auditory neural fibers 
to input of three tones. 



4.1 Model Parameters 

The parameters for ( |2.2|) are: cochlear length 3.5 cm; mass density 0.05 g/cm; 7 = 0.2; A = 
0.2 cm; N = 16.67 (dyne /cm 3 ); r$ = 0.001 g/(cm 2 ms). 

The middle ear serves as a low pass filter from sound pressure at eardrum to the displacement 
of stapes, and as a high pass filter from the sound pressure at eardrum to the acceleration. We 
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fit the data in [n| with the following gain factor 0^(0;) for each input e luJt : 

a M {u) = 30(4 + 0-0605^ 2 ((l - + ^w/wj 2 )- 1 / 2 ), (4-1) 
30 u 2 m 

where uj m = 4kHz, the middle ear characteristic frequency, and £ m = 0.7 the middle ear 
damping ratio. 

4.2 Dispersion Effect 

We illustrate the dispersion property of time dependent solutions numerically in case of one 
tone input. In top plot of Fig. [I], we show the profiles of computed traveling wave at t — 16 
ms with and without the selective damping term, e(x) ut- For the run, h = 0.01, k = 0.01, 
and x s = 3, j3 = 4. We see that for a tone input of 4kHz, 50dB, without selective damping, 
the BM displacement at t — 16 ms contains a long wave which manages to pass through the 
characteristic location in the interval (1.5, 2) and persists near the exit x = L. This is due to 
the dispersive long waves with slow decay as we discussed in the last section. In contrast, when 
the selective damping term is present, the long wave is damped as it moves close to the right 
end point x = L. As the selective damping term is only effective near x = L, the solutions in 
the interior of the domain are essentially the same. In the bottom plot of Fig. [TL we show the 
BM displacement at t = 32 ms, the dispersive long wave persists. In other runs (not shown 
here), we also observed growth in amplitude of long waves near x = L with time at a linear 
rate. The appearance of long waves seems to be independent of numerical methods, rather their 
existence is intrinsic to the model which is dispersive. 

4.3 Numerical Parameters 

For all our runs reported below, h — k — 0.01, e = 500, (3 = 4, x s = 3. The x s and (3 values 
should be properly increased (x s closer to x — 3.5) if there are low frequencies below 300Hz in 
the input frequencies. Smaller h and k values have been used to test that numerical solutions 
do not vary much with further refinement. The time step k should be made suitably smaller if 
input frequencies are as large as 16kHz for more resolution. 

4.4 Isodisplacement Curves 

For each input pressure wave of frequency ui, there is a unique location x^ G [0, L] where the 
maximum of BM displacement (peak in absolute value) is located. The i w is called characteristic 
place. Conversely, for each location x G [0, L] there is a frequency, called characteristic frequency 
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(CF) and denoted by CF(x), such that the BM response amplitude attains the maximum at 
x. In Fig. [I], we see that x^ G (1.5,2) for u = 4 kHz. The higher the amplitude of the input 
wave (denoted by A in ) , the higher the BM response peak. If the level of BM response is set at 
a fixed level, such as 1 nm, and a BM location x is given, we look for an input value Ai n so that 
|ii|(a;) = 1 nm. Here \u\ is the steady state BM response for input wave Ai n e lu)t . The plot of Ai n 
as a function of x gives the so called isodisplacement curve. The profile is asymmetric 

about the peak, and decays rapidly to zero beyond x^. So it takes higher A in to stimulate a 
point x > x^ to 1 nm than to the left. It takes least A in to stimulate j^Ka^) to 1 nm. So the 
isodisplacement curve has a minimum at x^, rises sharply at x > x u and gradually at x < x w . 

In Fig. |2|, we show model isodisplacement curves for five characteristic frequencies 500 Hz, 
2 kHz, 4 kHz, 6 kHz, 10 kHz. The threshold displacement is 1 nm. The curves showed the 
trend of sharper tips in high frequencies, wider tips for lower frequencies, and the asymmetry 
about the tips. The plot is comparable to the one in [23j. Isodisplacement curves are related 
to frequency selectivity of hearing. 



4.5 Two and Three Tone Interactions 

Two tone interaction is well-documented in the literature, here we illustrate that our model gives 
qualitatively the same results as shown in related figures in ||. We shall use BM displacement 
to detect the tonal interaction. Fig. ||| shows that the amplitude of a 5 kHz 70 dB tone drops 
(by nearly 50%) after the second tone of 2.4 kHz and 70dB is introduced, this reveals the so 
called low side suppression. 

In Fig. f|, we see that the 5 kHz tone at 50 dB (solid line) is suppressed (by about 30 %) after 
interacting with 6.7 kHz tone at 80 dB (plus-line), the so called high side suppression. The high 
side suppressor tone at 6.7 kHz has to be much higher in amplitude than a low side suppressor 
tone whose frequency is at least as far away from that of the suppressed tone. This shows the 
asymmetry of low and high side suppression as in || among others. 

In [[|, Fig. 7a showed cat auditory neural rate response at a fixed characteristic frequency 
(CF) demonstrating suppression effect for two tone (F 1; F 2 ) input, where F\ at fixed intensity 
takes three frequency values from below CF to above CF, for various values of the frequency 
F2 at 70 dB. In Fig. [|, we showed a qualitatively similar plot computed with our model, where 
CF=1 kHz, 50 dB Fi takes three frequency values (0.9, 1, 1.2) kHz, 70 dB F 2 increases its 
frequencies as 1.3 kHz, 1.4 kHz, 1.5 kHz, 1.6 kHz, 1.7 kHz, 1.8 kHz. 

The three tone interaction is much less documented in the literature, here we shall compare 
our model results qualitatively with the experimental data on auditory neural response to multi- 
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tones in [|J. Two tones are fixed at 4 kHz, 4.4 kHz, both at 50 dB. The third tone varies in 
frequency and in amplitude. In Fig. |6], we plot the BM response of two tones (4 kHz, plus 
line; 4.4 kHz, dotted line) as the third tone (line) increases its amplitude from 33 dB to 73 dB 
for three values of frequencies, 3.6 kHz, 3.8 kHz, 4.2 kHz. We see that the two tone response 
decreases for increasing amplitude of third tone; and that the suppression effect on the 4kHz 
tone is larger at 3.8 kHz (low side) than at 4.2 kHz (high side) for the same amplitude of third 
tone. The effect of third tone at 3.6 kHz is even larger than that at 4.2 kHz. In Fig. [7], we plot 
similarly the response of two tones as the third tone takes on frequencies 4.3 kHz, 4.6 kHz, 4.8 
kHz. The masking effect of third tone remains, except that the second tone at 4.4 kHz shows 
more difference as we observe the crossing of two tone curves when the frequency of third tone 
goes through that of the second tone (4.4 kHz). The third tone is more efficient in masking the 
second tone from the low side (4.3 kHz) than from the high side (4.6 kHz). The suppressing 
effect on the second tone is even weaker at 4.8 kHz. These features agree with those of Fig. 3 
of Deng and Geisler ||, e.g. compare frame 4 and frame 6 there. Strictly speaking the neural 
data in 0] are synchrony masking responses, and BM responses we computed correspond to 
rate masking responses. However, their data are relevant for qualitative comparison, as there is 
close (positive) correlation between the rate and phase behavior, see 0. 



5 Conclusions 

We discussed the dispersive property of the cochlear models for time domain computation, and 
slow decay of long waves. A selective longitudinal damping term has been introduced in the 
nonlinear nonlocal model studied previously jT^, [JJ. The new model is computed using a semi- 
implicit second order finite difference method in the time domain. We presented numerical 
results on isodisplacement curves, two tone suppression, and three tone interaction in qualita- 
tive agreement with earlier findings of Geisler || and experimental auditory neural data Deng 
and Geisler ||J]. These encouraging results prompt one to perform further study of masking 
mechanism in more complex tones and speech input based on our cochlear model. 
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Figure 1: BM displacement at t — 16 ms (top) and t = 32 ms (bottom), from a (4kHz, 50dB) 
input, with (line) and without (cross) selective damping. 
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Model Isodisplacement Curves with Threshold 1nm 
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Figure 2: Isodisplacement curves (1 nm BM displacement) at five characteristic frequencies 500 
Hz, 2 kHz, 4 kHz, 6 kHz, 10 kHz. 
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Figure 3: Low side masking effect of 2.4 kHz tone on 5kHz tone both at 70 dB. 



18 




0.4 



1 

r. 

/! 

i\l 


-- v« 

'i 

i 




- 5 kHz tone after interacting w. 6.7 kHz tone 


i i i 



0.5 1 1.5 2 2.5 3 3.5 

BM length in cm 

Figure 4: High side masking effect of 6.7 kHz tone at 80 dB on 5kHz tone at 50 dB. 
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Figure 5: The suppressed BM displacement (in nm) at characteristic frequency CF =1 kHz in 
the presence of two tone F 2 ) input, for 50 dB Fi = 0.9 kHz (dashed line), 1 kHz (circle-line), 
1.2 kHz (plus-line), as a function of F 2 which takes frequency values (1.3, 1.4, 1.5, 1.6, 1.7, 1.8) 
kHz at 70 dB. 
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Resp of two tones (4kHz,50dB),+— ;(4.4kHz,50dB),.— ; as 3rd tone ampl incrses 
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Figure 6: BM response of two tones (4 kHz, H — ; 4.4 kHz, dot-) as the third tone ( — ) increases 
its amplitude from 33 dB to 73 dB at three values of frequency, 3.6 kHz, 3.8 kHz, 4.2 kHz. 
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Resp of two tones (4kHz,50dB),+— ;(4.4kHz,50dB),.— ; as 3rd tone ampl incrses 
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Figure 7: BM response of two tones (4 kHz, H — ; 4.4 kHz, dot-) as the third tone ( — ) increases 
its amplitude from 33 dB to 73 dB at three values of frequency, 4.3 kHz, 4.6 kHz, 4.8 kHz. 
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